Including state-of-the-art physical understanding of thermal vacancies in Calphad models

A physically sound thermochemical model accounting for explicit thermal vacancies in elements and alloys is presented. The model transfers the latest theoretical understanding of vacancy formation into the Calphad formalism where it can extend currently available thermodynamic databases to cover vacancy concentrations without a complete re-assessment. The parametrization of the model is based on ab initio-calculated enthalpy of vacancy formation and two model parameters describing the excess heat capacity of vacancy formation. Excellent agreement is obtained with temperature-dependent vacancy concentrations and elemental heat capacities while reasonable extrapolation of phase stability to high temperatures is ensured. Extrapolation to multicomponent systems is reasonable and the long-standing Neumann–Kopp related problem in the Calphad community is solved since multicomponent solid solutions will no longer show fingerprints of elemental heat capacity peaks at their melting points. FCC-Ag, FCC-Al and FCC-Cu, FCC-Zn, FCC-Ni, BCC-Ti, and BCC-W are used as a demonstration, along with the Cu–Zn binary system.

www.nature.com/scientificreports/ formation on an atomic scale in two important directions: First, these studies challenged the previous consensus that the entropy of vacancy formation was temperature-independent 33 , which was based on the Arrhenius approximation, and showed that it is linearly dependent on temperature 33 . Second, different contributions to heat capacity have been quantified for several elements. While the required computational efforts prohibit the calculation of these contributions for multicomponent systems, for all elements studied high-temperature heat capacity changes (in first-order approximation) linearly with temperature if vacancies are ignored 30,[33][34][35] . Therefore, a Calphad model should be based on a temperature-dependent entropy of vacancy formation, and it is reasonable to assume that the heat capacity of the hypothetical vacancy-free metal changes linearly with temperature.
In this work, we combine the latest understanding of thermal vacancy formation from experimental and theoretical points of view into a thermochemical model that can be applied to and easily parametrized for multicomponent solid solutions up to very high temperatures. The thermochemical model is based on the compound energy formalism (CEF) and is fully compatible with existing CALPHAD databases but removes unphysical assignment of non-zero Gibbs energy of a hypothetical crystal consisting only of vacancies. It is also consistent with the recent conclusion that excess entropy of vacancy formation is temperature-dependent and zero at 0 K. The parametrization is reduced to two parameters, where one is accessible by computationally cheap calculations based on density functional theory and the other is fitted to experimentally determined vacancy concentrations. At the same time, it is ensured that the model extrapolates well to temperatures up to 6000 K, making it generally applicable to all solids. As a side effect, the long-standing Neumann-Kopp related problem in the CALPHAD community is solved so that multicomponent solid solutions will no longer show fingerprints of elemental heat capacity peaks at their respective melting points. As study systems, FCC-Ag, FCC-Al and FCC-Cu, FCC-Ni, BCC-Ti, and BCC-W are used. The model behavior for solid solutions is demonstrated for the FCC phase in the Cu-Zn system.
Due to its effects on materials properties, the formation of vacancies has been experimentally investigated for decades. Kraftmakher has compiled an excellent review of the known literature from an experimental point of view 6 . In pure metals, vacancies form because of the configurational entropy, while bond breaking leads to a positive vacancy formation enthalpy. In heat capacity measurements of pure metals, a non-linear increase in heat capacity is observed because with increasing temperature the vacancy concentration increases 6 . For a pure metal, Gibbs energy of vacancy formation G f (T) is considered the main quantity to describe vacancies and relates to the equilibrium vacancy concentration 7 , c eqm Va , by Eq. (1): Here, k B is the Boltzmann constant and T temperature. Calorimetric measurements are often relied upon to obtain c eqm Va values 6 . Despite the considerable impact vacancies have on mechanical and thermophysical properties, it is still rather difficult to accurately obtain accurate measurements of vacancy concentrations along the whole temperature range. Furthermore, experimentally measuring vacancy concentrations is a limited process that demands a complicated combination of different techniques 33 . The differential dilatometry (DD) method is considered the most accurate technique for measuring vacancy concentrations following the standards set by Hehenkamp 36 . Still, the results are most reliable between the melting point and about 80% of T melt only, as the results obtained below this point are too scattered to be reliable 6,33 . Equation (1) can easily be derived from a Gibbs energy minimization argument for a system consisting of one metal and vacancies, but it does not apply to solid solutions of metals or compounds such as (Ti, Al, Va)(N, Va) 5 and is completely inappropriate for phases such as (Ti, Va)(O, Va) and (V, Va)(O, Va) that exhibit vacancy concentrations on both metal and oxygen sublattice in the order of 20% 37 .
One of the downsides of experimental calorimetric measurements is that it is impossible to control vacancy concentrations. On the other hand, recent computational advances have enabled atomistic simulations that have significantly enhanced our understanding of vacancy thermodynamics 6 . Most importantly, Glensk et al. 38 and Gong et al. 31 have shown that the entropy of vacancy formation is temperature-dependent and reduces to zero at 0 K. This indicates that the Gibbs energy of vacancy formation must be described by a combination of enthalpy of vacancy formation and heat capacity of vacancy formation. Additionally, several ab-initio based studies have disentangled heat capacity contributions for vacancy-free metals 30,[33][34][35] . As is shown in Fig. 1, heat capacity of the vacancy-free metal can be approximated well by a linear temperature-dependence at high Temperatures (~ above 2/3 of the melting point) 30,[33][34][35] .
The CEF has been implemented on a wide scale in the Calphad method to describe thermodynamic properties of crystalline phases and to determine their phase equilibria 10,39 . With suitable parametrization, the CEF would be applicable also for solid solutions containing vacancies. In the last decades, several attempts have been made to model vacancies within the Calphad approach using the CEF, where it was necessary to assign a value for the Gibbs energy of a hypothetical pure vacancy solution endmember. While there is little debate that the Gibbs energy of pure vacancy (= nothing) should be equal to zero 7,19,40 , out of practical considerations Gibbs energy of pure vacancies is often described as non-zero, 60 kJ/mol 7 , 0.2 RT 17,20,21 , 10 RT 18 to 30 T 22 . The reason has been described by Franke 20 who provided a detailed mathematical analysis of the Gibbs energy of the solid solution for different parametrization, showing that an artificial local minimum in the Gibbs energy of the solution is generated for vacancy concentrations close to 100% due to the configurational entropy contribution if the Gibbs energy of the pure vacancy endmember is set to zero 20 . This artificial Gibbs energy minimum can be avoided if the pure vacancy endmember is sufficiently de-stabilized by assigning a positive Gibbs energy 20 . Another, rarely applied solution, is to limit the range of applicability of the solid solution description by setting a maximum vacancy concentration, e.g. 40% 22 . We will show below that both options -setting Gibbs energy of

Model description
In the CEF, Gibbs energy of an n-component solution containing vacancies is given by where y Va denotes the vacancy mole fraction, G Va is the Gibbs energy of a hypothetical pure-vacancies endmember, that is an empty crystal lattice occupied by vacancies only,y i and G i are the mole fraction and the Gibbs energy of the defect-free pure elements (1, . . . n) in the solution, respectively. R is the universal gas constant. L i,j E denotes the excess energy interaction parameter between two atoms i and j and L i,Va E between a host atom i and vacancy Va , both of which can be extended using Redlich-Kister (R-K) polynomials 41 . P denotes pressure and V represent the volume of the solution.
For simplicity, a pure metal A is considered here. The CEF extends the Gibbs energy of the solid solution between two end-members, a defect-free A crystal and a hypothetical vacancy-filled crystal 42 . For one mole of element A , the Gibbs energy is y A and G A are the mole fraction and the Gibbs energy of the defect-free crystal of A , respectively. The excess Gibbs energy parameter L A,Va E is where �H e (T) is the excess enthalpy and �S e (T) is the excess entropy. Note that both are assumed to be temperature-dependent.
Guan and Liu argued that it is necessary to consider the effect of pressure in the Gibbs energy model, to account for its impact on the total volume of the system, by including the PV term in a substitutional solution model that includes a vacancy 19 . Therefore, V is defined using the following formula: where L A,Va V is the excess volume interaction parameter between the host atom A and vacancy Va , extended using Redlich-Kister (R-K) polynomials 41 . Since we focus on Gibbs energy of vacancy containing metals and the PV term is small compared to the enthalpy of vacancy formation, we ignore second-order effects such as volume expansion and compressibility, for which more sophisticated models could however be used.
As introduced above, the choice for G Va has been under heavy debate. In many publications (e.g. Oates et al. 7 , Abe et al. 18 , Guan and Liu 19 , and Rogal et al. 8 , and) it is pointed out that a "crystal" without nuclei or electrons should collapse to zero volume for any pressure (> 0) and will contribute zero Gibbs energy and zero volume to a system. In that case, all thermodynamic properties of vacancy-containing crystals are described using the excess parameters ( . This is reasonable since the excess parameters describe the interaction between To tackle the mathematical problem presented by Franke 20 , i.e. the local minimum of the Gibbs energy for a vacancy concentration close to 100%, we consider both known options: First, we limit the range of applicability of the Gibbs energy model by the condition y Va < 0.2 , similar to a suggestion by Dinsdale et al. 22 . This limit is chosen to ensure that the very large vacancy concentrations observed experimentally, e.g. in Titanium and Vanadium monoxides in the range of ~2 0% on both sublattices of the Halite structure, are being covered. Note that it is only necessary to apply this limit for one sublattice (e.g. the metal sublattice of the Halite structure), while the vacancy concentration on the other sublattice (e.g. the oxygen sublattice of the Halite structure) may remain unconstrained (which enables the continuous description from Halite to FCC metal with metal vacancies and interstitially dissolved oxygen). Since most Calphad software packages do not offer such an option, we compare this to an arbitrarily chosen value of G Va = 2 J Kmol T . This contributes an unphysical entropy of vacancy formation at 0 K of −2 J Kmol which is compensated for by the excess Gibbs energy parameter L A,Va E which is still the most important parameter describing vacancy formation. This is possible since a crystal at 0 K exhibits an equilibrium vacancy concentration of zero which is perfectly in line with the ideal dilute limit.
Both parametrizations satisfy the condition that approaching the pure vacancy endmember from two different host crystals must result in the same lattice stability 43 and as we show in the supplementary Fig. S1 lead to equivalent descriptions for all vacancy concentrations of interest. All figures in the main part of this paper are identical for the two different model approaches.
Calphad databases rely on experimentally obtained heat capacity measurements to assess unary and higherorder systems. Since the effects of vacancy concentration cannot be suppressed during an experiment, the large enthalpy of formation of vacancies leads to a noticeable non-linear contribution in experimentally obtained heat capacity data. The G i term of Eq. (2) must however describe the defect-free crystal. This concern was raised by Abe et al. 18 , where the effects of vacancies on phase equilibria were accounted for twice. It was argued to be small and therefore ignored. However, it is clear that this is not the case for all elements 6 . Here, the implicit thermal vacancy effects are removed by adjusting the pure Gibbs energy description for the pure metal end members. According to Kraftmakher 6 , the thermal vacancies effects are limited under a certain temperature cutoff limit, usually 2/3T melt , where the their effects can be ignored. This is consistent with state-of-the-art atomistic simulations of high-temperature heat capacity of elements distinguishing harmonic, electronic and anharmonic contributions to heat capacity that indicate a rather linear high-temperature behavior, see e.g. W 34 , Al 34 , Ca 35 , and Ag 30 in Fig. 1. Therefore, the existing Gibbs energy descriptions 44 above T cutoff = 2/3T melt are modified to a linear extrapolation to exclude the implicit thermal vacancies effects.
In Calphad assessments, a constant value of the stable solid heat capacity description is usually adopted after T melt 10 . This avoids artificial re-stabilization of solid phases at very high temperatures due to too large heat capacities above the melting point 24,25,45 . Additionally, it is common in Calphad assessments to linearly interpolate heat capacity of multicomponent phases between the elements ("Neumann-Kopp rule"). This usually causes unrealistic heat capacity peaks for multicomponent phases at the melting points of the constituting elements 24,46 . To avoid the fingerprint from the elemental melting points in alloys' heat capacities we extend the linear extrapolation of elemental heat capacity beyond the melting point, however -to avoid re-stabilization at very high temperatures-only up to 2 * T melt , above which a constant value is adopted. Figure 2 shows the adopted heat capacity description for FCC-Al as example. The underestimation of heat capacity at the melting point is not a concern since at this temperature vacancies contribute to the experimentally measured heat capacity and will be described explicitly. where �H e (0K) is the excess enthalpy of formation at 0 K, and this quantity can be easily evaluated using ab initio calculations. �C e P (T) is the excess heat capacity of vacancy formation. The estimation of �H e (0K) quantity was carried out using multiple ab initio calculations for each investigated unary system (FCC-Al, FCC-Ag, FCC-Ni, BCC-W, FCC-Cu, FCC-Zn). The results were fitted using a simple quadratic function, yielding the value of �H e (0K) . This can be seen for FCC-Cu and FCC-Zn in Fig. 3. As can be seen here, all the obtained data points (up to 6.25% of vacancies) can be well described using a quadratic fit. Additionally, a comparison between the calculated �H e (0K) values from this work and corresponding experimentally determined values are shown in Table S1. The calculated values are in reasonable agreement with experimental ones.
�S e (T) is then defined as: where �S e (0K) is the excess entropy of formation at 0 K and is equal to zero 31,38 in case of G Va = 0 and to 2 J Kmol in the case of G Va = 2 J Kmol T which also leads to zero entropy of vacancy formation at 0K 31,38 . It was often assumed that vibrational entropy of defect formation was small and therefore ignorable 47 , which lines up with an Arrhenius-like description that presumes a constant entropy and enthalpy of formation. However, non-Arrhenius diffusion behavior was consistently observed in later experiments for several pure elements, such as Na 48 , K 49 , Ag 50 , and V 51 , and it was then determined that the Arrhenius approximation needed to be replaced by the Local Grüneisen Theory (LGT), which states that the temperature dependence of formation entropy can be approximated linearly for monovacancies and divacancies, which implies a temperature-dependent heat capacity of vacancy formation and scales in the Gibbs energy of vacancy formation by a T 2 approximation 33,38 . As will be demonstrated below, such a model extrapolates poorly to very high temperatures, so the following �C e P (T) expression was adopted: where a and b are fitting parameters. Setting b = 0 reduces the model to the LGT approximation but will result in constantly increasing entropy of vacancy formation and consequently when extrapolating the model to temperatures well above the melting point the solid phase will restabilize. In addition, it was recently shown by Grabowski et al. 26 that anharmonic vibrations and vacancies are strongly coupled and it can be expected that this coupling would increase further as vacancy concentrations increase in a hypothetical crystal superheated significantly above the melting point. Therefore, it is reasonable to assume that there is a non-linear �C e P (T) contribution at high-temperature and expand the function even further up to the second order in temperature. The choice of Eq. (8), makes it possible to re-write the �H e (T) and �S e (T) formulas as:  It should be noted that at low temperatures, the aT term dominates the entropy of vacancy formation while the bT 2 term is responsible for a well-behaving extrapolation to very high temperatures. The values of a and b are fitted to produce a final solution that satisfies two boundary conditions. (i) The first and hard condition is that the vacancy concentration at T melt matches the equilibrium vacancy concentration c eqm Va at T melt measured experimentally. This is achieved by fitting the a parameter in the first term of Eq. (8) with c eqm Va at T melt . (ii) The second and softer condition is that the resulting Gibbs energy description does not lead to re-stabilization at temperatures well above the liquidus temperature. This is achieved by determining the order of magnitude that the b parameter in the second term of Eq. (8) is required to satisfy this condition. To reduce the degrees of freedom during fitting, we suggest coupling the two parameters a and b . If b is set to zero unphysical re-stabilization of the solid phase happens at high temperatures, as shown in Fig. 4a applying the here-proposed methodology for FCC-Al while setting the value of the b parameter to zero. FCC-Al then re-stabilizes at around 3000 K. A model ignoring this will not be applicable to hard coatings of (Ti, Al, Va)(N,Va) 5 , which is in Calphad databases modelled as the same phase as the FCC phase and where at least the TiN endmember is stable above 3000 K.
Assigning a constant value to b applicable to all elements is not desirable since the −3bT 2 term in Eq. (10) should play a role only at the melting temperature, not at temperatures significantly below that. Therefore, we connect the b term to the a term by a coupling parameter α and the melting temperature, T melt using the following equation: So that the entropy of vacancy formation at the melting point is given by: Thus, the α parameter is a measure of how much the −3bT 2 contribution in Eq. (10) reduces the −2aT contribution to the entropy of vacancy concentration at the melting point. This procedure is in tradition of other works that have coupled Gibbs energy of vacancy formation parametrizations to the melting point 6 . As shown in Fig. 4a, for α = 0, re-stabilization of the solid phase occurs at high temperatures. For each element however, there is a lower limit of α that will ensure that the solid phase does not re-stabilize. Here, Al, Ag, Cu, Ni, Ti, and W have been studied to derive a universal coupling parameter α that ensures that artificial re-stabilization is avoided independent of the element. Using available equilibrium vacancy concentration c eqm Va data from the literature 6 , the value of the a parameter is fitted for different arbitrarily chosen values of α between 0 and 1 so that the value of the vacancy concentration at T melt would match the experimentally measured c eqm Va value without re-stabilization. All elements show re-stabilization below 6000 K for α = 0 as shown in Fig. 4a. For all FCCelements (Al, Ag, Cu, Ni) studied, the critical value of α is around 0.3, indicating that α FCC = 0.5 can be used as recommended universal coupling parameter for all FCC elements studied here. The BCC elements studied here, Ti and W, require a larger value, which leads us to recommend the value of α BCC = 0.75 adopted here. The reason for the larger value for Ti and W compared to the other elements is not clear and it is a task for future work to check whether this is a universally true for all BCC-metals. The interaction parameters for all assessed unary systems can be found in Table S2. Figure 4b shows that the Gibbs energy of FCC-Al extrapolates reasonable to www.nature.com/scientificreports/ 6000 K when the b parameter is set using the coupling parameter α FCC = 0.5 . No re-stabilization is observed. It must be noted that the number of fitting parameters is reduced to one by the universal setting of α . The single fitting parameter, a , is determined simply by the vacancy concentration at the melting point, c eqm Va . It is important to mention that the experimental c eqm Va data used were obtained using the Differential dilatometry (DD) method, which is regarded as the most reliable method of measuring equilibrium vacancy concentrations 6 .
The only parameter from Eqs. (2)-(13) not discussed extensively yet is the excess volume interaction parameter L A,Va V . As mentioned earlier, the importance of accounting for the excess volume was explored previously by Guan and Liu 19 , who accounted for the significant impact that the vacancy formation has on the total volume of the system. We use a somewhat simpler approach by accounting for the effect of vacancy concentration on volume in the interaction parameter L A,Va V . For this, V m was obtained from the same calculations used to determine �H e (0 K) . The excess volume is shown in Fig. 5 for FCC-Cu and FCC-Zn. In parameter studies we observed that the values obtained for V m contributed very little to the overall Gibbs energies of the examined systems and could therefore be ignored.
In the following paragraph, the model behavior is analysed and benchmarked against literature data.

Results and discussion
Here, the model performance is demonstrated with respect to heat capacity, vacancy concentration, Gibbs energy and entropy of vacancy formation, and the phase stability extrapolation to very high temperatures. Experimental data available in the literature that has not been used in the parametrization is shown.
Comparing with experimental C P data. For all the investigated unary systems, a comprehensive literature survey was carried out to collect as many heat capacity datasets as possible. These datasets will provide a visual guideline to ensure that the proposed description does not account for the thermal vacancies twice.
Having the highest melting point of all metals, BCC-W is considered a crucial case study for thermal vacancies. It represents the perfect example of how thermal vacancies can affect thermo-physical properties, such as heat capacity, in the vicinity of its melting point. For that reason, this element was investigated by Tang et al. 21 , who used the value, G Va m = 0.2RT , which was originally suggested by Franke 20 , who also chose to investigate the same BCC-W unary system in his work. It seemed ideal in this case to start with testing the proposed method on BCC-W and then compare the validity of the results accordingly. As shown in Fig. 6, the proposed model provided a very good description when compared with available experimental data for BCC-W and FCC-Ag. Furthermore, additional assessments were carried out for BCC-Ti, FCC-Al and FCC-Cu, and FCC-Ni, which are provided in the supplemental materials in Figs. S2-S5. Where available, the heat capacities of other descriptions including vacancies are included.

Vacancy fraction as a function of temperature.
Experimental vacancy concentration measurements obtained using the DD method for FCC-Al, FCC-Ag, and FCC-Cu are compared against the results produced by the proposed model in Fig. 7. The vacancy concentration at the melting point was used during the parametrization, so perfect agreement at this temperature is not a surprise. However, the model also matches the experimental data not used during parametrization accurately, especially close to T melt , where relative experimental errors become smaller due to the larger magnitude of vacancy concentrations, making measurements easier. This provides further support for the model's validity and physical integrity.
Gibbs energy and entropy of vacancy formation. The vacancy concentration data from Fig. 7 can be converted to Gibbs energy of vacancy formation, which was used by Glensk et al. 38 to benchmark their atomistic simulations against the experiments.     30 . Entropy of vacancy formation is shown in Fig. 9 for FCC-Ni and BCC-W. For Ni, the model reproduces the trend obtained from atomistic simulations 31 qualitatively, by reproducing a positive temperature dependence of entropy of vacancy formation at low temperatures and a reasonable order of magnitude at the melting point. For W, temperature-dependent entropy of vacancy formation is not available. However, the comparison to existing Calphad models 20,21 clearly shows the different approach suggested here: The unphysical assignment of negative entropy of vacancy formation is avoided in the model suggested here.
High-temperature extrapolation. All remaining unary systems were assessed using the proposed model. The results show that the proposed model was capable of maintaining sensible phase stability in all cases and can be found in supplementary materials, Figs. S6-S11.
Avoiding Neumann-Kopp artifacts in solid solutions. The model's heat capacity description managed to avoid Neumann-Kopp artifacts in solid solutions when compared with the existing SGTE description. Heat capacity for Cu-20% Zn is shown in Fig. S12.
Binary and higher-order systems extrapolation. The FCC phase in the Cu-Zn system was extended to the model proposed here while keeping all Cu-Zn interactions the same. The resulting phase diagram does not change significantly, as shown in Fig. S13.  www.nature.com/scientificreports/

Methods
The estimation of �H e (0K) and V m was carried out using multiple ab initio calculations for supercell structures containing 15, 31, and 47 atoms for the FCC metals and 15, 23, 35 and 53 atoms for BCC-W, respectively, plus one vacancy. The calculations were performed through the VASP software using GGA pseudo-potentials. Volume and atomic positions have been relaxed. Thus, �H mix (0K) was calculated for each investigated unary system (FCC-Al, FCC-Ag, FCC-Ni, BCC-W, FCC-Cu, FCC-Zn). From the same calculations, V m was obtained accordingly.

Conclusion
The proposed thermochemical model joins the most recent experimental vacancy formation findings with the latest state-of-the-art physical understanding of vacancies and their effects. It is shown that the long-standing debate in the Calphad community on the choice of G Va can be settled since both options ( G Va = 0 and G Va = positive but small) can result in equivalent descriptions if parametrized using the methodology proposed here. Other Calphad-related concerns, such as extrapolation to multicomponent systems and re-stabilization of solids at high temperatures, were addressed successfully in this work. Additionally, the model results shows excellent agreement with experimental heat capacity and vacancy concentration data as well as theoretical entropy of vacancy formation data. The proposed model can easily be applied to other elements or compounds if experimental and ab initio data are available. Explicitly accounting for thermal vacancies contributions in multicomponent solutions is thus possible and can be very beneficial when modelling processes like diffusion or precipitation.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.